【题解】 [ZJOI2014]力 FFT bzoj3527/luogu3338

这个式子有点……乱。

嗯，我们来推一推式子……推一推式子。

原式推一推，那么就是：

E_{i} = \frac{F_{i}}{q_{i}} = \sum_{j = 1}^{i} \frac{q_{j}}{(i - j)^{2}} - \sum_{j = i + 1}^{n} \frac{q_{j}}{(i - j)^{2}} ​

$E_i = \frac{F_i}{q_i} = \sum_{j=1}^{i} \frac{q_j}{(i-j)^2} - \sum_{j=i+1}^{n} \frac{q_j}{(i-j)^2}$

令 $x = \frac{1}{y^2}$ ，那么：

E_{i} = \frac{F_{i}}{q_{i}} = \sum_{j = 1}^{i} q_{j} x_{i - j} - \sum_{j = i + 1}^{n} q_{j} x_{j - i} ​

$E_i = \frac{F_i}{q_i} = \sum_{j=1}^{i} q_j x_{i-j} - \sum_{j=i+1}^{n} q_j x_{j-i}$

还可以写成：

E_{i} = \sum_{j = 1}^{i} q_{j} x_{i - j} - \sum_{j = i + 1}^{n} q_{j} x_{j - i} ​

$E_i = \sum_{j=1}^{i} q_j x_{i-j} - \sum_{j=i+1}^{n} q_j x_{j-i}$

令 $S_i = q_{n-i+1} $ ，那么式子变成了：

E_{i} = \sum_{j = 1}^{i} q_{j} x_{i - j} - \sum_{j = i + 1}^{n} p_{n - j + 1} x_{j - i}

$E_i = \sum_{j=1}^{i} q_j x_{i-j} - \sum_{j=i+1}^{n} p_{n-j+1} x_{j-i}$

这个时候我们可以发现， $\sum_{j=1}^{i} q_j x_{i-j}$ 和 $\sum_{j=i+1}^{n} p_{n-j} x_{j-i}$ 都是卷积，那么我们可以跑两遍 $FFT$ ，分别求出上面的两个式子，记录为 $A,B$ 。最后的答案就是 $A[i].x - B[n+1-i].x$ 了。

FFT不用做太多修改，套模板跑就行(本来就是模板)。

Code:

#include<cstdio> 
#include<cmath>
#include<string>
#define ll long long
#define inf 0x3f3f3f3f
#define RI register int 
#define PI 3.1415926535898
using namespace std;
const int N=4e5+2;
int n,limit=1,filp[N];
template <typename _Tp> inline _Tp max(const _Tp&x,const _Tp&y){return x>y?x:y;} 
template <typename _Tp> inline _Tp min(const _Tp&x,const _Tp&y){return x<y?x:y;}
template <typename _Tp> inline void IN(_Tp&x){
    char ch;bool flag=0;x=0;
    while(ch=getchar(),!isdigit(ch))if(ch=='-')flag=1;
    while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
    if(flag)x=-x;
}
struct complex{complex(double a=0,double b=0){x=a,y=b;}double x,y;};
complex operator + (complex a,complex b){return complex(a.x+b.x,a.y+b.y);};
complex operator - (complex a,complex b){return complex(a.x-b.x,a.y-b.y);};
complex operator * (complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);};
complex A[N],B[N],C[N];
inline void FFT(complex *f,short inv){
    for(RI i=0;i<limit;++i)if(i<filp[i]){complex t=f[i];f[i]=f[filp[i]],f[filp[i]]=t;}
    for(RI p=2;p<=limit;p<<=1){
        RI len=p/2;
        complex tmp=complex(cos(PI/len),inv*sin(PI/len));
        for(RI k=0;k<limit;k+=p){
            complex buf=complex(1,0);
            for(RI l=k;l<k+len;++l){
                complex t=buf*f[l+len];
                f[l+len]=f[l]-t,f[l]=f[l]+t,buf=buf*tmp;
            }
        }
    }
}
int main(){
    scanf("%d",&n);int cnt=0;
	for(RI i=1;i<=n;++i){scanf("%lf",&A[i].x),B[n+1-i].x=A[i].x;}
	for(RI i=1;i<=n;++i)C[i].x=(1.0/double(i))/double(i);
	while(limit<=(n<<1))limit<<=1,cnt++;
	for(RI i=0;i<limit;++i)filp[i]=((filp[i>>1]>>1)|((i&1)<<(cnt-1)));
	FFT(A,1);FFT(B,1);FFT(C,1);
	for(RI i=0;i<=limit;++i)A[i]=A[i]*C[i],B[i]=B[i]*C[i];
	FFT(A,-1);FFT(B,-1);
	for(RI i=0;i<=limit;++i)A[i].x/=limit,B[i].x/=limit;
	for(RI i=1;i<=n;++i)printf("%.3lf\n",-B[n+1-i].x+A[i].x);
    return 0;
}

本文标题:【题解】 [ZJOI2014]力 FFT bzoj3527/luogu3338

文章作者:Qiuly

发布时间:2019年02月14日 - 00:00

最后更新:2019年03月29日 - 13:54

原始链接:http://qiulyblog.github.io/2019/02/14/[题解]luoguP3338/

许可协议: 署名-非商业性使用-禁止演绎 4.0 国际转载请保留原文链接及作者。